;; Non-robust linear algebra functions
;; Copyright (C) 2021  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.linalg.tridiagonal)

(defn- back-substitution
  "Perform the back substition of a tridiagonal system of linear equations.
  It is assumed that the input sequences are already reversed."
  [b c d]
  (loop [xs (list (/ (first d) (first b))) ; Use a list to accumulate in reverse
         b (rest b)
         c c
         d (rest d)]
    (if (seq d)
      (let [x (/ (- (first d) (* (first c) (first xs))) (first b))]
        (recur (cons x xs) (rest b) (rest c) (rest d)))
      xs)))

;;; Tridiagonal matrix algorithm:
;;; https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
(defn tridiagonal
  "Solve a tridiagonal system of linear equations."
  [a b c d]
  (loop [b' (list (first b))     ; Use a list to accumulate in reverse
         d' (list (first d))     ; Ditto
         a a
         b (rest b)
         c' c
         d (rest d)]
    (if (seq d)
      (let [w (/ (first a) (first b'))]
        (recur (cons (- (first b) (* w (first c'))) b')
               (cons (- (first d) (* w (first d'))) d')
               (rest a)
               (rest b)
               (rest c')
               (rest d)))
      (back-substitution b' (reverse c) d'))))
